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ABSTRACT 

Numerical models of the dynamically extended atmospheres of long period vari- 
able or Mira stars have shown that their winds have a very simple, power law structure 
when averaged over the pulsation cycle. This structure is stable and robust despite 
the pulsational wave disturbances, and appears to be strongly self-regulated. Observa- 
tional studies support these conclusions. The numerical models also show that dust- 
free winds are nearly adiabatic, with little heating or cooling. However, the classical, 
steady, adiabatic wind solution to the hydrodynamic equations fails to account for an 
extensive region of nearly constant outflow velocity. An important process or group 
of processes is missing from this solution. Since gas parcels moving out in the wind 
are periodically overrun by pulsational waves, we investigate analytic solutions which 
include the effects of wave pressure, heating, and the resulting entropy changes. 

In the case of dust-free winds we find that only a modest amount of wave pressure 
is needed to derive an analytic model for a steady, constant velocity, locally adiabatic 
outflow. Wave pressure is represented with a term like that in the Reynolds turbulence 
equation for the mean velocity. The waves damp relatively quickly with radius, as 
a result of the work they do in driving the mean flow. Although the pressure from 
individual waves is modest, the waves are likely the primary agent of the self- regulation 
of the dust-free winds. 

In dusty Miras, the numerical models show the radiation pressure on grains and 
the subsequent momentum transfer to the gas, play the dominant roles in driving 
the wind, and wave pressure is not very important. In the models of the dusty wind, 
the gas variables also adopt a power law dependence on radius. Heating is required 
at all radii to maintain this flow, and grain heating and heat transfer to the gas are 
significant. Both hydrodynamic and gas/grain thermal feedbacks can transform the 
flow towards particular self-regulated forms. 
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1 INTRODUCTION 

The astrophysical theory of hydrodyna mic wind s bega n with 
the solar wind eq uations proposed by IParkerl y.958) , also 
see <IParkerill963l) . Parker's model predicted that the pres- 
sure gradient between the sun's outer atmosphere and the 
interstellar medium would drive the wind. In subsequent 
years Parker's solutions to the hydrodynamic equations have 
proven to be a powerf ul tool for stu dying stellar winds in 
general. (The analogous iBondl ( 1952) solutions for spherical 
accretion problems have proven to be equally important.) 
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iBirdl (1964a be) constructed shock heated models of 
the solar corona and wind, which can be regarded as 
Parker winds with an additional heating term. However, 
Bird's models were not designed to capture the highly 
nonlinear shock dynamics of long period variable star at- 
mospheres (henceforth LPVs). The large-amplitude, quasi- 
ball istic motions beh ind LP V shocks were not considered 
fsee lWillson fc Hilll979l and lHill fc WillsoiJl97^ . Nonethe- 
less, a number of important results relevant to Miras were 
anticipated, including: acoustic heating of the circumstellar 
gas, a dynamic balance between shock heating and expan- 
sion cooling in the wind, and self-regulation in the wind flow. 
Bird proposed these as driving forces for the solar wind and 
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corona, which are now believed to be the resu lt of magneto - 
hydrodynamic dissipation (e.g., the review of|Parkerj[l997). 
At the present time, these processes appear more relevant 
to the warm extended at mospheres and winds of dust-free 
LPVs (see the reviews of lWillson et~al]ll997l Iwillsonll200d 
for discussions of these structures) . As will be described be- 
low, self-regulating processes are very important, but Bird's 
conjecture that the self-regulation works to maintain a con- 
stant Mach number throughout the flow is not strictly cor- 
rect. 

Hart mann and MacGregor j Hartman n fc MacGregorl 
1980, Hartmann & MacGregor 1982, also see 



iMacGreeor fc Charbonneaul Il997l) proposed an appli- 
cation of solar corona and wind models with Alfven wave 
heating and pressure to LPV winds. These processes 
are similar to the shock pressure and heating processes 
described below. The Hartmann and MacGregor models 
assumed an approximately hydrostatic corona, with waves 
viewed as a modest perturbation, in contrast to the strong 
shocks observed in LPVs. Specifically, these are polytropic 
wind models, without shock entropy generation. Both shock 
and Alfven waves may play a role in LPV winds, but we 
believe that the strong pulsational shocks are dominant, so 
we will not consider magnetic effects in this paper. 

Because LPVs have very large amplitude pulsa- 
tions, they might be expected to provide difficulties for 
smooth wind solutions of the hy drodyn a mic e q uation s. In- 
deed, the numeri c al mo dels of iBowera <ll988l . Q99CT). and 
iBowen fc Willsonl (I1991T) show that the atmospheres of 
pulsating Miras are highly dynamic. This is also true 
of all pu blished numerical mo d els of these and relate d 
stars, e.g.jFleischer et all Jl992l).lFeuchtinger et alJ <1993l). 
Hofner fc Dorfil lll997lL ISteffen. Sz czerba 1 fc Schoe nbemed 
Jl998ft . lWinters et al] | |200ClL and lHofner et alH2003l) . Most 
of these models do not include non-LTE radiative cooling, 
so in the remainder of this paper we will primarily refer to 
Bo wen's models, and th e particular runs described below. 
See IWillson et al] 1119971) for a comparison of various mod- 
els.) 

In the numerical models, the density near the photo- 
sphere retains a roughly exponential decline, but departs 
from this to a power law dependence at density value that 
depends on the mass outflow rate. In the transition from 
exponential to power law decline, the dynamical lifting by 
shocks followed by (in some cases) radiation pressure on 
grains lifts and accelerates the material. In this flatter den- 
sity profile, more gas has been lifted to large radii, and 
this along with radiation pressure on dust grains in some 
cases, provide the means for driving much enhanced winds. 
Bowen's models show that the inner atmosphere is a complex 
region where shock accel eration and heating and cooling all 
play important roles f see IWillson et al. 1997J). On the other 
hand, the wind region at large radii has a simpler structure, 
and the gas variables averaged over pulsation period can be 
well approximated as simple power law functions of radius 
(sec Figures 1-3 and discussion below). Such simple profile 
forms suggest that it should be possible to construct simple 
analytic models in the tradition of Parker and Bird, at least 
for the smooth wind region. 

We will not consider LPV observations in detail in this 
paper, since our main goal is to understand the regular- 
ities revealed by numerical models. In addition, the con- 



straints on theory provided by observation (e.g., spectra) 
are indirect, and generally the comparison between theory 
and observation is best don e with the aid of detailed nu- 
merical models (see e.g.. IWillson 120001 . iTei et al . 2003 and 
references therein). Nonetheless, some recent observations 
provide quite direct information about the gas variables in 
the winds and extended atmospheres of LPVs, and a brief 
mention of these provides a context for the subsequent work. 

First of all, a major assumption of both analytic 
and numerical models is that the LPV winds are approx- 
imately spherically symmetric. In spite of evidence for 
modest deviations from symmetry, or asymmetries in bi- 
nary systems, there is much evidence from a variety of 
wavebands for approximat e symmetry in most cases (e.g., 
iBuiarrabal fc Alcoleairi99lT) . 

Recent infrared interferometric observations of "dust 
shells" around Miras are also generally consistent with 
spherically symmetric models, though not with uniform dust 
density distributions. The data are better fit by models in- 
cluding a few distinct she l ls where the dust emission is high 
(see e.g..lHale et al.lll997llLopez et al]ll997L iMonnier et alJ 
1997, Fong, Meixncr, & Shah 2003, and references therein). 



However, these observations probe a region located at a 
distance of only a few stellar radii, where conditions are 
described as dynamic and complex. In fact, the numeri- 
cal models show that this is the region where shocks have 
grown to very nonlinear amplitudes, where dust forms if it 
is able to, and where the wind is just beginning to be ac- 
ce lerated (the "shock acceleration" and "dissipation" zones 
of IWillson et al.lll997l) . Thus, these observations generally 
support the picture provided by the models, and do not 
contradict the notion of a generally smooth wind structure. 
Observations and models are in general agree- 
ment about basic wind parameters like the mass loss 
rate or flow velocity. From CO observations (e.g., 



I Buiarrabal fc Alcoleal Il99ll iKahane fc Jural Il99l lYoune 
1995, iKerschba um fc Olofssonl Il998l iG roenew egen et al. 
Il999l and lWinters et al.ll200l also see lAlard et al.ll200ll for 
mid-infrared results) we know that mass loss rates range 
from 10 -7 Mo/yr up t o 10~ 5 Mp,/yr for Miras in agreement 
wi th Bowen's mo dels <BowerJll98a IBowen fc Willsonlll99ll. 
fc IWirlsonll2000D . Above about W~ 5 M Q /yr the stars are 
generally classified as OH-IR sources. In this paper we are 
most interested in low mass loss rates, in relatively dust-free 
cases. 

The dust free models develop an approximately con- 
stant outflow velocity with a speed well below the escape 
speed, and this pattern persists over many stellar radii. In 
this region the flow is also subsonic. Thus, this constant ve- 
locity is not the same as the coasting flow in standard wind 
models well beyond the sonic point. The observations also 
find low wind velocities (of order 5-10 km/s) and little ev- 
idence for changes in the flow velocity through the wind. 
For the Miras CO observations are able to probe the flow 
at large distances from the star, and so, provide important 
constraints on any model. Interestingly, Cepheid wind ve- 
lociti es also seem to be low relative to the escape speed 
( Sass elov fc Les ter 1994aj|bJ). 

In the following sections we will derive analytic wind 
models, described by hydrodynamic equations which includ- 
ing approximate, averaged terms, like those in the Reynolds 
turbulence equations, to represent the effects of nonlinear 
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Figure 1. Density profile from a dust-free, hydrodynamical model like those of lBowerJ il988f) . The model was run for a very long time 
to allow a large steady wind region to develop and relax out to r ~ 10 15 cm. (See text for details). The dotted line shows an a p oc r~ 2 
function (normalized to the numerical model at log(r) = 14.5) for comparison. 
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Figure 2. Radial velocity profile from the end of the same dust-free Bowen model shown in the previous figure. Comparison to the 
dotted line at zero velocity shows the low, constant value of the outflow in the steady wind region. The upper dashed line shows the 
local escape velocity. The lower dashed curve shows a fit derived from an analytic model, see Sec. 2.6 for details. 
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Figure 3. Temperature profile from the end of the same long, dust-free Bowen model shown in the previous figures. The dotted line 
shows an a T oc r _1 function (normalized to the numerical model at log(r) = 14.5) for comparison. 
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pulsational waves on the flow. We then discuss the self- 
regulatory feedbacks in these winds, and use several ap- 
proaches to argue that these processes select constant ve- 
locity outflows over other possible wind solutions. 



2 WAVE PRESSURE IN DUST FREE WINDS 

2.1 Classical Steady Wind Solutions Compared to 
a Numerical Model 

The usual Parker wind equations are derived from the 
steady, spherically symmetric hydrodynamic continuity and 
momentum equations, which can be written as, 



M — Airr pu 
du 1 dP 



GM 



dr 



p dr 



(1) 

(2) 



where the variables, M, M, p, P, u, and r denote the mass 
loss rate, total mass of the star, the density, pressure, wind 
speed, and the distance from the center of the star, respec- 
tively. Classical steady wind solutions also assume a poly- 
tropic equation of state, P = Kp 1 , with constant K, 7. In 
this paper we will assume a constant rate of mass loss as 
well. 

Two properties of the outer region of the atmosphere 
suggest that it can be described as approximately locally 
adiabatic: 1) the temperatures there are in a range where 
radiative cooling is inefficient, and 2) the shocks in this re- 
gion are relatively weak. The latter feature is evident in 
Figures 1-3, which show density, velocity, and temperature 
profiles taken from the end of a very long run (about 1000 
pulsation cycles) of Bowen's modeling code, without dust. 

Detailed descriptions of the atmospheric dynamics can 
be found in Bowcn ( 1988). Here we merely note that the code 
numerically solves the one-dimensiona l, spherically symmet- 
ric, hy drodynamic equations with a iRichtmver fc Mortonl 
( 1967) type Lagrangian method, including a conventional 
artificial viscosity algorithm for stabilizing strong shocks. 
A number of different chemical rate, momentum and en- 
ergy sources are calculated explicitly for each gas element. 
These include: the atomic radiative cooling rates, the free 
electron abundance, energy changes due to ionization and 
recombination. They also include simple approximations to 
the radiative transfer, grain formation, radiation pressure on 
grains, and grain heating and cooling. The effects of most of 
these processes will not be considered in this paper, which fo- 
cusses instead on understanding the global hydrodynamics. 
However, the numerical models provide a standard for com- 
parison, as well as assurance that the analytic models are not 
unrealistic. The parameters of the numerical model shown 
in Figures 1-3 are essentially the same as those of a dust- 
free, I.OMq, "standard" model of Bowen ( 1988), though the 
model shown here was computed with increased spatial res- 
olution. The numerical models will be compared to analytic 
models below, especially in Section 2.6. 

The wind density profile shown in Figure 1 has an ap- 
proximately 1/r 2 form, while the form of the corresponding 
temperature profile is nearly 1/r (Fig. 3). Note that beyond 
a radius of 10 cm in the figures the model wind flow is 
not fully relaxed. (Indeed, Fig. 3 suggests that the flow is 
not thermally relaxed beyond log(r) = 14.5.) There is, in 



fact, a polytropic solution to equations (1) and (2), with a 
1/r 2 density profile, and a 1/r temperature profile. However, 
the value of the polytropic index in this model is 7 = 3/2. 
Not only does this not correspond to the adiabatic value of 
7 = 5/3 expected in the absence of cooling and heating, it 
does not have any obvious physical meaning. 

We should note, however, that Parker argued against 
values of 7 > 3/2 on the grounds that larger values did not 
yield "a solution beginning at low velocity close to the sun 
and e xtending outward to zero pressures at infinity" JParkerl 
1963, pg. 61). He also noted the need for a heating source 
when 7 < 3/2. 

What about the adiabatic, 7 = 5/3 wind solution? To 
illustrate the problems with this classical solution we inte- 
grate equation (2) from an inner radius n to an arbitrary 
outer radius r, use the ideal gas law and employ the defi- 
nition of the adiabatic sound speed to derive the following 
Bernoulli equation, 



Ml 



7" 



2 2 

c — ci 



Mel 



(3) 



The numerical models show that the sound speed c decreases 
with radius approximately as r~ ' , as does the escape ve- 
locity M e . Thus, according to the equation (3), the flow ve- 
locity u must decrease with radius by comparable amounts. 
It does not; Bowen's models show constant velocity in the 
wind at all times (see Fig. 2). 

We note also that, if all the wind flow lies on a sin- 
gle adiabat in the thermodynamic phase space, and has a 
1/r 2 density profile, then the ideal gas law implies that the 
temperature goes as T oc r -4 ' 3 . With the long run we have 
performed, which yields a large wind region, the numerical 
profile does appear shallower than an adiabatic one. 

Maintaining a shallower temperature profile requires ex- 
tra heating and momentum sources in the wind. In fact, 
a likely reason for discrepancies between numerical models 
and the adiabatic solution is that, although the pulsational 
shocks are weak in the wind region, their momentum, heat, 
and entropy inputs are not negligible. 



2.2 Not Quite Adiabatic Solutions with Wave 
Pressure and Heating 

2.2.1 Wave Pressure 

The question is how to capture the effects of waves in a rel- 
atively simple analytic description of the overall wind flow, 
averaged over many pulsation cycles? The problem is like 
that of describing turbulent flows in which large fluctua- 
tions on many scales coexist with a mean flow and long-lived 
coherent structures. The characteristics of the regular, pul- 
sational shocks, which propagate outward through the wind, 
are very different from the stochastic, fluctuations described 
in th e theory of well-developed, homogeneous turbulence 
(e.g.. iTennekes fc Lumlevlll972l lMcComblll990l) . Nonethe- 
less, in both cases we are dealing with disturbances on small 
spatial scales and short timescales, which we don't necessar- 
ily need to resolve in order to describe their average effect 
on the flow. 

In fact, the fundamental idea of turbulence theory, that 
the flow can be divided into a mean and a fluctuating part, 
provides a very appropriate approach to generalizing the 
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classical stellar wind equations. We cannot simply adopt the 
averaging procedures used to derive the Reynolds equations 
of turbulence theory, since these depend on the random na- 
ture of the fluctuations in homogeneous turbulence. We can, 
however, provide physical justifications for including similar 
'averaged' fluctuation terms in a set of phenomenological 
Reynolds equations for mean wind flows. 

The numerical models show that each pulsational shock 
provides a compression and an outward boost to an indi- 
vidual gas element, and the gas element f ollows a quasi- 
ballistic free-fall trajectory behind the sh ock ( Wi llson fc Hill 
| l97flllHill fc Willsonlll979i lBowedll98sL and the review of 
IWillson et alJll997f) . More precisely, for a periodic pattern 
of motion in the atmosphere, the material at larger r must 
have smaller outward postshock speeds. This picture is good 
deep in the atmosphere. If the postshock speed is slightly 
too high, the mass element encounters the next shock at 
slightly larger than it did the last one. This gives a small 
net outward motion, as is seen in the atmosphere around 
1.5-2.5 stellar radii. Where the flow dominates, the material 
still responds nearly ballistically, setting up a stable pattern 
relative to the mean flow speed and providing a small in- 
crement of momentum to the gas in each cycle. This allows 
us to separate the mean flow from the oscillations, which 
average approximately to zero (see Appendix A for details) . 

To incorporate the effects of the shocks on the mean 
flow, an additional term is included in the (inviscid) mo- 
mentum equation. With the approximation of the previous 
paragraph, this term can be written as the divergence of a 
Reynolds stress, which in the spherically symmetric case is 
just the radial derivative of the mean square velocity fluctu- 
ation, or d(a 2 )/dr. Physically, this is the (one-dimensional) 
wave pressure gradient. In the Bernoulli equation (3) it con- 
tributes a a 2 term, like the sound speed term. 



2.2.2 Shock Heating 

The net effect of shock heating can be included as a term in 
the mean, steady state energy equation (see Appendix A), 
which can be written as, 



1 _ 2 8 I , :; 

2" d~r^ rpU 



dr 
GMpu 
r2 



r pue 



+ Vp. 



+ r~ 



dr 



(4) 



The variables, e and Y denote the internal energy per 
unit gas mass, and the shock heating function, respectively. 
Henceforth, we define u as the sum of the mean (U =< u >) 
and fluctuating (a — (u — U)) velocity components. In tur- 
bulence theory, we assume that this is an ensemble average 
over many nearly identical systems, distinguished by the de- 
tails of their random fluctuations. In the present case the 
average is also assumed to be over all pulsational phases. 
Similarly, p, e, and P generally have both mean and fluctu- 
ating parts, but we will not need to adopt special symbols 
for the separate parts (except in Appendix A). 

We re-emphasize that in this equation we neglect the 
effects of the interaction between wind material and the stel- 
lar radiation, the interaction with the radiation generated in 
different parts of the winds, and energy exchanges between 
gas and dust grains. 



The Reynolds formalism provides an additional equa- 
tion for the mean square fluctuating velocity (see McComb, 
sec. 1.3.2). In the case of a steady, inviscid, constant mean 
velocity radial flow, this equation can be written, 



or 



d_ 
dr 

-2a' 



< (u-Uf >+- < (u-U)P> 



dU_ 
dr 



■2I\ 



(•») 



The first two terms on the right-hand-side represent turbu- 
lent (or wave) energy diffusion by nonlinear couplings, which 
we expect to fall off quickly with radius in this case, and so 
will neglect them. (This neglect is part of a closure approx- 
imation for the Reynolds moment equation set.) What re- 
mains is a relation between the local wave pressure gradient 
and the shock heating, 



2 dr 



(•") 
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+ C7 



,dV_ 

dr ' 



(6) 



(where the last term is small in winds with nearly constant 
mean velocity). 

Although we now have the required energy and momen- 
tum source terms and a relation between them, the equation 
set is not quite complete. 



2.2.3 Equation of State 

We employ a locally, but not globally, adiabatic equation of 
state to describe the effects of shock heating in the outer 
atmosphere. Physically, in each pulsation cycle, any par- 
cel of gas goes through an irreversible (non-closed) cycle 
in a p-V (pressure - specific volume) phase diagram. A 
pulsational shock wave pushes a gas parcel along a Hugo- 
niot curve in the phase diagr am, and off its initial adiabat 
(e.g. Landau & Lifshitz 1959). Radiative cooling at constant 
pressure behind the shock then moves it to lower specific vol- 
ume (i.e., compresses it). Downstream from the shock the 
parcel evolves along a new adiabat to higher specific vol- 
umes and lower pressures in the rarefaction region, until it 
is hit by the next shock. A possible laboratory analogue of 
this physics is provided by the reverberating shock cavities 
produced in the (gun-type) experiments (see lHolmes et alJ 
Il995l) . 

For weak shocks the specific entropy changes are not 
great, but they accumulate. Gas parcels farther out in the 
wind have experienced more of these cycles, and so, lie on 
adiabats that are farther from the initial one than those in 
the inner wind. We assume that the specific entropy of the 
mean flow is time-independent and a smooth function of ra- 
dius, like the other gas variables. We believe this gradual 
entropy change is an important factor in making the wind 
region of Bowen's dust-free Mira models appear quite adi- 
abatic, but with a constant outflow velocity, which is not 
characteristic of a classical adiabatic wind. 

We can write the locally adiabatic equation of state as, 



K{r)p^ 



(7) 



where 7 = 5/3. Because gas parcels at different radii are on 
different ad iabats, K is a function of the r adius, rather than 
a constant. IlKlahr fc Bodenhe imcr (2003J) have also recently 
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considered model Keplerian disks with entropy gradients due 
to local heating.) 

Since the gas is locally adiabatic, the perfect gas rela- 
tionship for e is valid (locally) at any radius, 



(8) 



P 7- 



This relation completes the equation set. We have in- 
troduced three new terms into the Parker wind equations, 
which describe wave pressure, shock heating, and shock en- 
tropy production. None of these extra phenomenological 
terms are found in the gas equa tions that ar e the basis of 
published numerical models (e.g.. lBowenll988f) . Those mod- 
els follow the time-dependent pulsation and shock phenom- 
ena explicitly. We believe all of the phenomenological terms 
are necessary to describe the locally adiabatic mean flow in 
pulsationally driven winds. 

We note that the classical, linear treatment of waves 
propagating into a stellar atmosphere indicates that long 
period waves should be reflected at the surface. Pijpers and 
collaborators, described the effect s of stochastic ac oustic 
waves in detail in a series of paper s iPijpers fc Hearn | l989a. 
[ Pijpers fc Hearnl 1989 j. iKoninx fc Piipersl I1993J IPiiperj 
Il993l . and iPiipere Il995l) . They found that such waves can 



jropagate in t he co ntext of a general outflow (see iPiipersI 
1993).|BowerJ al99CJ) found that the power needed to main- 
tain large amplitude photospheric oscillations is actually 
lower for long periods, and these induce nonlinear wave prop- 
agation (shocks) traveling out through the atmosphere. Here 
we take the existence of transmitted waves as a given, and 
focus on the large-scale, nonadiabatic effects of mildly non- 
linear waves in the wind region only. 



2.3 Constant Velocity Winds and Entropy 
Production 

Motivated by the observations, we begin by considering the 
simple case in which the mean wind velocity is constant, 
U = Uq. In this case, equation (1) immediately gives p, and 
the mean flow momentum (eq. (2)) can be written, 



ldP 

p dr 



GM d a 

—5 — -j-(<7 =0, 

H dr 



(9) 



where p and P are now the Reynolds average quantities. 

With the additional assumption that all three terms on 
the left-hand-side of this equation scale the same way with 
radius throughout the wind flow, we can readily obtain a 
solution to this and the energy equation. This assumption 
appears quite strong, but is physically reasonable for the 
following reason. The wind extends over a very large range 
of radii. The gas thermal energy and temperature do not fall 
off as rapidly with radius in dust-free models as expected 
in an adiabatic flow. This suggests a heat source operates 
throughout; shock heating is the only available source. (We 
note that photo-heating of grains will have a similar scaling 
with radius, but we assuming negligible grain populations 
for the present.) This fact ties the thermal pressure to the 
shocks and suggests that it scales like the wave pressure. 

Beyond this physical argument, we can briefly note the 
mathematical possibilities. First is the possibility that ther- 
mal pressure dominates in part of the wind and wave pres- 
sure elsewhere. In the former region we will have a Parker 



thermal wind, with the observational difficulties mentioned 
above. In the wave dominated region, we will get the same 
scaling with radius, which derives from balancing the gravi- 
tational term. Since shock compression does provide thermal 
heating, this limit isn't physically self-consistent. A second 
possibility is that the thermal and wave pressures have dif- 
ferent functional forms whose sum contrives to balance the 
gravitational term. For example, they could both have the 
same power law form derived below, but modulated by os- 
cillatory parts which have the same amplitude for both, but 
which are 180° out of phase. Any other solution of this type 
would have to be similarly fine-tuned. Given that the wind 
extends over orders of magnitude in radial extent, this seems 
physically contrived. We will argue later that the solution 
below is also preferred by self-regulatory processes. 

Specifically, we assume that a 2 = A/r for some constant 
A, then equation (9) becomes, 



dP 



- (gm - a\ ■ 



-dr. 



Now p can be substituted from the continuity equation 
and the result can be integrated (from a radius r to 00) to 
get 



(10) 



p= l(GM-A) 

3 r 

assuming p and P go to zero at infinite radius. It is interest- 
ing that this equation has the form P oc p 3 ' 2 of a constant 
velocity Parker wind, since i oc p', though we here view 
the solution as only locally adiabatic, rather than globally 
barotropic. 

Next, we want to simplify the energy equation (4) in 
this constant velocity case. We need the following result in 
taking the mean of the term on the left-hand-side of equation 
(4), 



pu > = < p 



U 3 +3a 2 U 



(11) 



where p is used here as the total density (mean plus fluctu- 
ating parts). This equation assumes that the mean of odd 
powers of the fluctuating velocity vanish, and that the cor- 
relation < p(u — U) > is negligible (see Appendix A). Hence- 
forth, we drop the brackets <>, and assume again that p 
and P refer to the mean quantities. 

When we use this result, pull Uq factors out of partial 
derivatives, and substitute for e from equation (8), the en- 
ergy equation becomes, 



3 / pUo \ d 
2 I r 2 I dr' 



:(rV) + |(y) 

'7-2^ Uo d V P ) GMpUo 



dr 



7 — 1 J r 2 dr 



+ Tp, 



(12) 



To further simplify we can substitute for P with equa- 
tion (10), for A with equation (6), use the substitution 
a 2 = A/r. The result reduces to, 



A = —GM, 



and yields in turn, 

5 GM 

^22^ — ■ 



(13) 
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The physical interpretation of the last equation is that 
the shock heating rate is doing work against the gravita- 
tional potential GM/r at each radius, on a flow timescale 
of about r/Uo- It is reducing the effective gravity (i.e., the 
factor GM-A in eq. (10)), that must be overcome by flow 
down the thermal pressure gradient. 

Combining equation (10) with equation (7) and the con- 
tinuity equation gives the following expression for the vari- 
able K in equation (7). 



K{r) = -GM\ 



M 
4tt[/ 



and that equation becomes, 



P=^GM I" 
11 \ 4tyU 



r-3 p' 



(14) 



(15) 



a generalized (position-dependent) polytropic relation. 

To better appreciate why K varies with radius, recall 
that the entropy (per unit mass) of a perfect gas is, 



S = log(P/p 



(16) 



so equation (7) implies S = log(K). Thus, K is directly re- 
lated to entropy, and its radial dependence implies a mean 
radial entropy gradient. As any gas element flows outward, 
it is overtaken by shocks moving through the wind, and each 
shock increases the entropy of the element. At any fixed 
radius the entropy production is balanced by the outward 
transport. 

The value for the constant A above implies that the 
wave energy is about 45% of the gravitational potential en- 
ergy at all radii. This implies that a should be about half 
the value of the escape velocity at any radius. The velocity 
"jitter" due to shocks in Bowen's dust-free numerical mod- 
els is considerably less than this, generally less than a third 
of the escape velocity. 

However, it is not clear that this velocity jitter should be 
directly identified with a, see section 2.6 below. We expect 
difficulties in capturing the full range of velocity variation in 
the models. Most of this variation occurs across thin shocks, 
and these short wavelength variations are not well resolved 
numerically. 

2.4 Accelerating Winds 

The constant velocity wind described above is an especially 
simple solution to the equations. In this section we consider 
steady accelerating winds with radially dependent mean ve- 
locities. We restrict our consideration to power law similar- 
ity solutions. Mathematically, we are essentially exploring a 
specific kind of variation of the previous solution, but self- 
similar winds are the most relevant physically. 
Specifically, we assume U = Uo(r/ro) 

™ r „-(2+i) 



47riin 

is, 



so p = 
In this case, the mean momentum equation 



dU _ _ldP_ 
dr p dr 



GM 



d 
1U- 



(17) 



To solve this equation we will again assume a specific 
form for the wave pressure term, that is, 



A 



B 



H) 



(18) 



where the first term is as in the constant velocity case, and 
the second term is designed to adjust the wave pressure in 
accordance with the flow acceleration. Like A, the coefficient 
B is independent of radius, but does depend on the exponent 
8. Note that in the case 8 — 0, this expression does not 
coincide with the previous one for constant flow unless we 
add an additional term of -B. This term complicates the 
algebra of the model considerably, and since we usually have 
\B\ « A/ro (see below), we will neglect it for the moment. 

We could also add other power law terms to the expres- 
sion for the wave pressure, or pursue a series solution for the 
wind variables. We believe this would only complicate the 
equations without adding much physical insight. 

An expression for the pressure can now be derived by 
integrating the momentum equation as before. We obtain, 



P(r) = 




GM-A 



\ p(r) 






+ 



i +%)*Hk 



(19) 



The function K(r) can be derived by substituting equation 
(7) for P. 

Next we substitute equation (18) (with nonconstant U) 
into the a 2 equation (5) (including the term for the flow of 
energy from the mean field to the fluctuating velocity field, 
— 2a ^=r). This produces a rather complex expression for the 
heating, 



r = 



2r„ 



V J r \r \ r\) 



35-1- 



(20) 



Finally, the above expressions can be substituted into 
the energy equation to derive the following expressions for 
the coefficients A and B as functions of the exponents 7 and 
8 by equating the coefficients of the equal powers of 8, 



(2- 7 )-(3 + <5)(7-l) 

A = - — ^ J — -GM, 



(2-7)-(2-5)(3 + 5)( 7 -l) 



B = 



"-(7- 


-l)(2-<5) + (2- 7 )2<5 


5(7- 


-1)(2- 8) -45(2-7) 



u$. 



(21) 



(22) 



The equations (20)- (22) specify the net heating required 
to maintain the power law velocity profile. They do not di- 
rectly take into account the dependence of specific heat- 
ing processes, like shock heating, on the gas quantities, but 
rather assume that such processes can be regulated to the 
form above. However, the expression for the heating profile 
(eq. (20)) is so complex that it seems unlikely that realistic 
shock heating and radiative cooling processes would gener- 
ate it. This is in contrast to the simple heating profile of 
equation (13) for the constant velocity outflow, and sug- 
gests that the simpler form would be preferred (as it is in 
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Bowen's numerical models). This question will be examined 
more carefully in the next section. 

Before continuing, however, we note a couple of pecu- 
liarities of the equations above. First of all, in the limit of 
\5\ « 1 and 7 = 5/3 (locally), the expressions for A and B 
become, 



A -rM 1 + i 6 



B 



l^-To 5 



(23) 



(24) 



B is negative, and not directly proportional to 5. This means 
that in the accelerating wind, the wave pressure could be 
reduced at r < ro relative to the constant velocity flow. 
Physically, it is more reasonable to simply identify ro as the 
inner radius of the wind (see eq. (18)). 

Another peculiarity is that if, as we expect, the heating 
is dominated by the first term in equation (20), it becomes 
negative when 5 is increased above a value of 1/2, which 
implies a cooling process. This also seems unphysical. 



2.5 Self- Regulation 

2.5.1 General Considerations 

The simple appearance of the profiles of gas variables in 
Bowen's dust-free numerical models suggests that the sys- 
tem relaxes to a self-similar solution, independent of the de- 
tails of the shock heating, and subsequent expansion cooling. 
We have suggested above that one key to deriving this solu- 
tion is the assumption that the gas in the wind is only locally 
adiabatic. This assumption allows for entropy production in 
weak shocks, which gradually evolves a gas element through 
a sequence of adiabats as it moves outward in the wind. 

The thermal state of the gas determines the environ- 
ment traversed by the shocks, which in turn are responsible 
for producing this state. This feedback in the shock heat- 
ing mechanism eventually works to regulate the wind to the 
1/r 2 density profile. 

The process of self-regulation is well illustrated by the 
early relaxation that occurs in Bowen's numerical models, 
which are started with small (but increasing) pulsation am- 
plitudes, and with an exponential atmospheric density pro- 
file. As described in IBowenl 119881 the traveling wave parts 
of the pulsations steepen into strong shocks as they propa- 
gate down the initial exponential atmospheric profile. These 
shocks launch gas pa rcel s out on quasi-ballist ic trajectories 
l|Willson fc Hilll979l and lHill fc WiHsonll979h to such large 
radii that they are not able to fall back to their starting 
points before being hit by the next shock. As a result, the ini- 
tial profile is stretched out into power law form. This causes 
a feedback, the shock acceleration is decreased in the flatter 
pressure profile, and the ballistic launch velocities are mod- 
erated. Thus, if the pressure profile is steeper than 1/r 3 , 
then the shock amplitude grows, pushing material out, until 
it matches that form (see Appendix B). If the density gradi- 
ent is shallower, however, the loss of energy in pushing the 
gas will lead to declining shock amplitudes, and it will not 
be possible to hold the shallow profile up against gravity. A 
pressure gradient that is neither too steep, nor too flat, also 
maintains the flat outflow velocity. 



The steady gas variable profiles that are eventually es- 
tablished with finite amplitude pulsations include another 
region between the exponential atmosphere where shocks 
grow nonlinearly, and the power law wind region. This is the 
shock dissipation region where the power law density profile 
flattens, but the temperature and pressure increase with ra- 
dius to a maximum, forming a warm "Caloris phere" and the 
inner boundary of the wind. See the review of lWillson et alJ 
( 1997) and references therein for more details. The develop- 
ment of this region in the numerical models provides a clear 
example of the effects of shock dissipation when density and 
pressure profiles flatten. 

This qualitative discussion on relaxatio n processes, 
which was partially anticipated by 'Birdl ll964ol . helps us un- 
derstand why there is a preferred wind profile. However, it is 
not yet sufficiently precise to predict the form of that profile. 
Next we will consider some more quantitative approaches to 
the problem. 



2.5.2 Nonequilibrium Thermodynamics and Least 
Dissipation 

The relaxation to the simple, constant velocity flow is 
driven, in part, by the tendency of such a system towards 
a state of least dissipation, or minimum entropy produc- 
tion. Such states often correspond to those of maximum 
entropy content. The Theorem of Minimum Entropy Pro- 
duction for steady, nonequilibrium states derives from work 
of Helmholtz, though it has been generalized and used in 
a variety of applications in the last few decades (see e.| 
iGlansdorff fc Prigogindlllml |Prigogindll98cl . IWoodslll99e 
Proofs of the different versions of this theorem involve sub- 
stantial restrictions, e.g., that the system is not too far from 
thermodynamic equilibrium, or that there are linear rela- 
tions between the flow variables and their driving force or 
sou rce terms (like the familiar linear stress-strain relation), 
see lWoodslJl99d) . 

These theorems have not been used much in astro- 
physics, where there are many non-equilibrium systems, but 
where these systems are often highly time-dependent. On 
the other hand, the winds of long-period variables appear 
to provide an astrophysical situation where the theorem can 
be usefully applied. If pulsations begin with small ampli- 
tudes and build up steadily, then the wind can adjust to 
this changing driving, and never find itself too far from a 
set of local equilibria. Moreover, equation (9) (neglecting 
the last term) appears to give the linear relat ion required o f 
the dissipative quantities in some proofs (see|Woodsj[l996). 
The minimization of shock dissipation also seems to be in 
accord with the qualitative considerations of the previous 
subsection. 

We can apply the theorem by comparing the net pro- 
duction of the entropy in the accelerating wind models. The 
theorem implies that the wind with minimal entropy pro- 
duction is the preferred state. In principle, we need to con- 
sider the details of the dissipational processes, which is very 
complex in the present case. Fortunately, it seems clear phys- 
ically that the entropy production in the steady winds must 
be directly related to the shock heating function pr (see eq. 
(20)). Thus, for the purposes of making a simple estimate, 
we define a function, H, as follows, 
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47ir pTdr 



(25) 



where ro and n are adopted inner and outer bounds of the 
wind region. The coefficient (2/M) is included merely to 
absorb that common factor from the continuity equation 
into the definition of H, which then has units of velocity 
squared. Substitution for pT and integration yields, 



H 




+ 2B(S) 1 



2« 



(26) 



where the functions A and B are given by equations (21) 
and (22) with 7 = 5/3, and xi = ri/ro- Note that the outer 
boundary of the wind, n, will generally be a function of 6, 
which governs the density falloff. For example, we can make 
a simple estimate by assuming that the outer boundary is 
where the mean flow velocity U equals the escape velocity. 
In that case, we have, 



(27) 



The dimensionless factor (U 2 ro)/GM, which appears 
both in the ratio of coefficients in equation (26) Bro/A and 
in equation (27), is the primary similarity parameter of the 
problem. For convenience, we name it, 




Z = 



Upro 
GM' 



(28) 



Now, using equations (21), (22), and (27) we can ap- 
proximate H as, 



H 



GM (1 




(29) 



Then we take the derivative with respect to <5, and set 
it equal to zero to find the dissipation extrema. The result- 
ing expression is very complicated, and not worth recording 
here without some simplifying manipulations and approxi- 
mations. 

As yet, the only extra approximation we have used in 
deriving the expression for H (or its derivative) is the outer 
boundary estimate of equation (27). However, at this point, 
it is helpful to adopt the approximation that Uq << 2M ; 
i.e., that the inner-edge wind velocity is much less than the 
inner-edge escape velocity. This assumption is clearly satis- 
fied in the numerical models. Since B cxUn, and — oc QM. 
this is equivalent to the assumption that |B| << A/ro al- 
ready noted above. However, it is not wise to neglect the 
B-term in equation (26) above if 8 is small and negative, 
because the term x\ can be relatively large. 

Next, consider limits on the magnitude of 8. When 8 
is positive, the first term in equation (29) is larger than the 
second, and for 8 > 1/2 it is negative. Then the dissipation H 
would be negative, which is unphysical. Thus, we do not need 
to consider large, positive values of <5. Similarly, when -^ < 
8 < -^ the first term makes H negative, and unphysical. For 



large negative values of 8 the second term of equation (29) is 

negative and dominant. In sum, it seems that the physically 

relevant region is where the exponent \S\ < 1/2. 

These approximations justify dropping terms of order 
1 
Zq or Zq 1 + 26 , and then the derivative of H is given by, 



-^ = = (m 2 - 140(5 + 100 H 8<5 3 + 128 2 + (v) - I ) \ 



# 



-8<5 4 - 4<5 3 + 128 2 - 1945 + 32 



+ 



2 2<T - 28 + 11 49<T - 140(5 + 100 x 



(-«"-" + «) f 



-2.S 

T+2J 



32(5 



('-')( 



2S 2 -2<5 + ll x 



(-w +M ) f 



26 

I 2:S 



(30) 



This equation can be solved as a quadratic in the vari- 

25 

able (Zo/2) 1 + 26 as a function of 8. The value of Zq itself is 
then obtained for each value of 8. We find that equation (30) 
yields no positive, real values of Zq for positive values of 8 
less than about 1/3. With the constraints cited above this 
means that there is no significant range of positive values of 
8 that give physical solutions. That is, solutions with wind 
velocity increasing with radius do not satisfy the least dissi- 
pation constraint. On the other hand, there are solutions to 
equation (30) for negative values of 8. The value of Z in- 
creases monotonically (and ri/ro decreases monotonically) 
as 8 goes from to increasing negative values. If we make 
the reasonable requirement that the radial extent of the out- 
flow with u < v e , that is ri/ro, be greater than a few, then 
\8\ must be less than about 0.1. As the magnitude of \8\ is 
decreased below 0.1, ri/ro rises very rapidly, so the wind is 
very large for such small values of 8. 

In sum, the Least Dissipation Theorem, as applied to 
the shock heating function, has constrained the family of 
wind solutions. Specifically, the parameter 8 is limited to a 
small range for optimal winds. Because the physically rele- 
vant values of 8 have a very small magnitude all the mem- 
bers of this one-parameter family are nearly constant ve- 
locity winds. This result is not completely rigorous, though 
physically reasonable. It is worth emphasizing that because 
self-regulation is a global process, global constraints deter- 
mine the wind structure. 



2.5.3 Constraints from the Bernoulli Equation 

As an integral of the momentum equation, the Bernoulli 
equation provides another global constraint on wind struc- 
ture. For the accelerating winds with wave pressure de- 
scribed equation (17), the corresponding Bernoulli equation 



UHr 1 + lir)_ +(j2 ,1 

2 7' - 1 y ' 2 



Cs, 



(31) 
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where, as in equation (3), v e is the escape velocity, 7' is the 
global effective ratio of specific heats, and the integration 
constant Cs is a function of the wind acceleration exponent 
5 only. 

In the case of an accelerating wind we have expressions 
for all of the terms on the left-hand-side of this equation, 
except the sound speed term. Using equation (9) or (19) 
above, we find, 



c 

7 



7 



GM- 



(2 + 6)r 



•2D 



U$r 



GM - A 



1+2,5 



Since for 
3 + <5 



3 + <S 

small, we have P oc p^+s , 



and then, 



(32) 



(33) 



U , the first term on the right-hand-side of this equation 
will generally dominate the others. (Note that in the limit 
of small |<5|, the last term is also proportional to 5, so it does 
not exceed the first term.) However, the sign of the first term 
is such that it would imply a falling U for a positive S, and 
a rising U for negative 8, which contradicts the definition 
U oc r . The only alternative is 8 — 0. 

As before, this conclusion is not mathematically rigor- 
ous, because it is not always true that the first term on 
the right-hand-side of the above equation must dominate. 
For example, the parameter B goes to positive infinity as 
S approaches a value of 10/7. However, such cases are ex- 
ceptional, and unphysical. In sum, we generally expect that 
thermal winds driven by nonlinear acoustic waves have con- 
stant outflow velocities. 



The value of the constant Cg can be determined by 
evaluating the left-hand-side of the Bernoulli equation at 
the inner edge. In Bowen's numerical models of dust-free 
winds we find that, 



eg 



r 



1 



2c « 



2 



(34) 



The models also indicate that uo ^ Co, and that Uo is much 
less than either of these. Thus, the terms on the left-hand- 
side of equation (31) nearly cancel, and \C's\ has a relatively 
small value. 

Next, we compare an accelerating wind, with \8\ small, 
but nonzero, to the constant velocity wind having the same 
outflow velocity at the inner wind radius. We also assume 
that the inner radius and stellar mass are the same for both 
winds. Taking the difference between the Bernoulli equations 
for the two winds we obtain, 



i^-i*)-. p^ Vf- 



.^w; 



2 + 6 "''■■" ' 




B = Cs-C s = . (35) 



To determine the difference between constants on the 
right-hand-side of this equation we evaluate it at r — ro- 
We substitute the result back into equation (35), and solve 
this equation for the fractional difference in the mean flow 
velocity, as a function of r, 



U 2 - ul 




(36) 



Because Uq is generally much less than the escape ve- 
locity throughout the wind, except possibly the outermost 
parts, and because the magnitude of \B\ is comparable to 



2.6 Wave Pressure and the Outer Wind 

One remaining aspect of the dust-free winds to consider is 
what happens in the outermost regions? For brevity, we will 
confine our discussion to the constant velocity wind. The 
escape velocity v e , the sound speed c (eq.(32)), and the ve- 
locity dispersion associated with the waves, a (eq. (18)) all 

-1 
decrease with as r 2 in this case. In the outer wind the con- 
stant flow velocity will eventually equal and exceed each of 
these velocities in turn. 

The effective a is generally somewhat less than the 
sound speed c, so we might expect a transition from the 
constant velocity wind to a nearly adiabatic outflow before 
the material escapes. Specifically, for the constant velocity 
wind with 7' = 3/2, the parameter A = 2GM/5 (eq. (21)), 
and equations (18) and (32) yield, 



a ~ 0.45i> e o 

C ~ 0A7VeO- 



(37) 



Thus, we see that when Uo equals cr(r), it very nearly equals 
c. However, these velocities are less than half the escape 
velocity, and there is no reason to expect an abrupt increase 
of c(r). That is, we do not expect a transition to supersonic 
flow through a Parker critical point (where Uo = c = v e ). 
Physically, as Uo becomes comparable to a and c, and all 
about equal to v e /2, we might expect successive pulsational 
waves to propel material above the local escape velocity, and 
thus free of the star. It appears that the outermost wind 
must be nonsteady. 

The result that the sound speed is about half the lo- 
cal escape velocity is consistent with the numerical models. 
However, the velocity amplitude of the pulsations is much 
less than that. As discussed above, this likely due in part to 
the limited spatial resolution of the models. Moreover, the 
wave pressure used here could well represent the aggregate 
effect of several pulsational waves. (In fact, it is possible for 
one shock to overtake and merge with another, though this 
process is not well-resolved in the numerical models.) 

We can illustrate this last point with a more detailed 
comparison to the numerical model. We consider what wave 
pressure boosts are sufficient to push a gas element through 
a radius change of order unity, Sr/r = 1. The wave pres- 
sure drives the constant mean flow, and the time for the 
gas element to flow this distance is 8t = Sr/U. The num- 
ber of pulsational shocks passing through the gas element 
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in this time is about, TV = St/T, where T is the pulsation 
period. We identify the net wave pressure with the sum of 
the pressures in each shock to get, a = NAv . According 
to equation (37) the model predicts, a ~ 0.45i> e , and so it 
also predicts, 

Av ~ ^LeO- (38) 



N 

As a specific example we take r = Sr — 10 14 cm, and 
U — 0.6 km/ s. This is a radius where the numerical model 
shows that strong shocks present at smaller radii have set- 
tled into the steady wind. There are still significant velocity 
fluctuations well beyond this radius, but the adopted mean 
outflow velocity is a fair average for the whole wind region. 
With these values and a pulsation period of about one year, 
we compute the following values for the quantities defined 
in the previous paragraph, 



St = 1.67 x 10 J s, N ~ 54, v e 
and Av ~ 1.1 km/s. 



18. km/s, 



(39) 



We set Av equal to this latter value at the adopted radius, 
use the scaling of the constant velocity outflow model [p oc 
r -1 ' 2 ), and add the mean outflow U to derive the lower 
dashed curve in Figure 2. This curve is generally comparable 
to the velocity jumps of the shocks in the inner wind, where 
these jumps are fairly well resolved. In fact, the velocity 
jumps appear to be somewhat larger than the curve in the 
inner wind, but this is reasonable since some of the shock 
energy will go into thermal heating as well as doing work 
against the star's gravity. 

In the present example the work per mass of gas moved 
is —GM/(2r). The corresponding shock heating per gas 
mass, integrated along the path, is F5t — 0.11 (GM/r). 
The integrated work done by the wave pressure is about 
half what is needed. This confirms the impression given by 
the figure that there is not too much wave energy available 
for heating. In sum, in terms of the wind driving energetics, 
wave pressure and thermal pressure play comparable roles. 
A fraction of the thermal energy is replaced by shock dissi- 
pation, as required for an effective adiabatic index of 3/2. 



3 SELF-REGULATION IN RADIATIVE AND 
DUST DRIVEN WINDS 

If dust grains, atoms or molecules are able to absorb a frac- 
tion of the substantial Mira radiative output, then radiation 
pressure can play a substantial role in driving or accelerating 
th e wind. In the case of typical Mira variables , the models 
of lBowenl (Il98gl) . and iBowen fc Willsonl Jl99ll) . show that 
models with significant radiation pressure on grains have 
higher mass loss rates and higher wind velocities. The wind 
regions are also much smaller in these models because the 
wind velocity surpasses the escape velocity at a smaller radii. 
The pulsational hydrodynamics is essentially frozen out 
in a dusty wind. There is, however, another regulatory pro- 
cess that limits the wind velocity and keeps the flow from 
over-refrigerating. In fact, the numerical models show that 
the gas temperature tracks the radiative equilibrium tem- 
perature quite closely, which implies a coupling of photons, 
grains and gas atoms. Atom-grain collisions a re the mecha - 
nism of both momentum and energy transfers (|BowerJl988). 



The gas flow is driven by the momentum transfer from each 
grain to individual atoms, and these are inevitably accompa- 
nied by energy transfers. The interaction will also decelerate 
the grain, and will also result in grain ablation if the rela- 
tive velocity is too great. Ablation will reduce the grain cross 
section, and the radiative driving. Slow wind speeds allow 
more time for grain growth in dense regions, leading to more 
radiative driving. 

In particular, once the grains relax to a fixed size dis- 
tribution, the fraction of the photons they absorb, and the 
momentum input from them, decrease as 1/Y 2 . In this case, 
two terms dominate the right-hand-side of a momentum 
equation, like equation (2), the gravitational term and the 
radiative driving source term. These terms scale the same 
way with radius, allowing a constant velocity outflow solu- 
tion. Other solutions to the momentum equation are possi- 
ble. However, because of the regulatory feedbacks, we expect 
that the process of grain growth towards an equilibrium size 
distribution should naturally be accompanied by wind ac- 
celeration to a constant terminal velocity. 



4 CONCLUSIONS 

In summary, the numerical models indicate that the winds 
that develop in the extended atmosphere of long-period vari- 
ables have a simple power law structure. The isentropic 
Parker wind solution that matches this structure has a 
barotropic index of 7 = 3/2. In a real gas heating is required 
to maintain this index, but the Parker solution does not pro- 
vide information about the heating source. In the case of 
non-dusty Miras, where radiation pressure is unimportant, 
we believe that the power law structure is generated and 
supported by shock waves which travel through the wind. 
Since shocks dissipate energy and generate entropy, these 
winds have significant heat inputs and entropy gradients. In 
section 2 we presented phenomenological equations includ- 
ing the three relevant terms: wave pressure, wave heating 
and large-scale entropy gradients. We then studied a family 
of analytic thermal wind solutions to these equations with 
power law velocity profiles. These wind solutions generalize 
the classical isentropic Parker solutions to cases in which the 
gas is only locally adiabatic. The simplest member of this 
family has a constant outflow velocity, and matches numer- 
ical models of non-dusty Mira winds quite well. 

Section 2 concluded with a discussion of why the con- 
stant outflow velocity solution may be preferred in nondusty 
Miras. Both the Least Dissipation Theorem of linearized, 
nonequilibrium thermodynamics, and the Bernoulli equa- 
tion with reasonable physical constraints indicate that this 
is a preferred solution among the family of power law winds. 
This result provides reassurance of the basic correctness of 
Bowen's numerical models, and the validity of mass loss es- 
timates predicted by them. 

In Section 3 we considered dusty winds in long-period 
variables. In this case the winds are driven by radiation 
pressure on dust grains, and the flow "freezes out" in the 
sense that thermal, wave-driven or turbulent pressures are 
not dynamically important. This allows the wind velocity 
to attain supersonic speeds without passing throught the 
critical point of classical thermal wind theory. Nonetheless, 
Bowen's numerical models show that there are significant 
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couplings between the radiation field, an equilibrium distri- 
bution of dust grains, and the atomic gas. These couplings 
allow a regulated, constant velocity outflow to form and be 
maintained. 

We believe that the model of dust-free Mira winds may 
provide a paradigm, and that gas dynamic solutions with 
acoustic wave pressure may have much more general appli- 
cation in astrophysics. Specifically, such solutions will be 
relevant whenever there is substantial supersonic or magne- 
tosonic (velocity) noise on small scales as compared to large 
scale velocities and thermal pressure gradients. Laboratory 
examples ar e provided by t he rev erber ating shock caviti es 
described in iHolmes et"al] Jl995l) and IWeir et"aH i 1996(1 1. 
Other possible astrophysical examples include: 1) accretion 
disks with strong acoustic-convective turbulence, 2) galac- 
tic winds or outflows driven by continuing supernova shocks, 
and 3) the hot gas halos of galaxy clusters, containing galax- 
ies moving at transonic or supersonic velocities, and into 
which there is continuing infall. The defining property is 
that the entropy or the function K have a nonzero gradient. 
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APPENDIX A: APPENDIX A: DENSITY AND 
VELOCITY AVERAGES 

In this appendix we give a brief, qualitative discussion of how 
the density and velocity fluctuations driven by pulsational 
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shocks propagating through an LPV wind can be divided 
into a mean and a fluctuating part that averages to nearly 
zero over a pulsational period. We also discuss the derivation 
of the Reynolds equations for this flow. To begin, consider 
the trajectories of a couple of adjacent gas elements relative 
to the local mean flow. We choose the gas elements that are 
sufficiently close together that the local mean outflow can be 
approximated as constant velocity, even if it is not globally 
constant. 

Next, we divide the pulsational cycle into two parts. The 
first part begins when the gas element is hit by a shock. This 
impact boosts the fiducial element outward, and after a brief 
delay does the same to the neighboring element at slightly 
larger radius. They are pushed closer together in the shock 
compression, so both op and <5u > 0, i.e., greater than the 
values of the mean flow (see e.g., Fig. 2 of lBowenlll988l . 

The fiducial gas element is closer to the star than its 
neighbor, and so, has a shorter free-fall time. As a result of 
this and the fact that it received its outward velocity burst 
first, it will decelerate ahead of its neighbor. This ballistic 
description is equivalent to the statement that the shock is 
followed by a rarefaction wave. We define the beginning of 
the second part of the cycle as beginning at the moment 
when the velocity is reduced to the mean flow value, and 
continuing until the next shock impact. The density is re- 
duced to the mean value at about the same time. Thus, 
during this second interval we have, 5p, and Su < 0. 

Because the pulsational shocks are not very strong in 
the wind region, the magnitudes of 5p, and Su are not large, 
and the cycle is nearly symmetric in the sense that the two 
intervals are about equal to half a cycle. Thus, the net den- 
sity and velocity fluctuations over one cycle are small, i.e., of 
higher order than St/t, the fractional radius change in that 
time. This justifies the Reynolds procedure of writing each 
gas quantity as a sum of a mean and a fluctuating part. E.g., 
for the velocity u — U + a, and for the other fluid variables, 
p =< p > + p', P =< P > + P', e =< e > + e, etc., where 
the primes denote the fluctuating parts. 

The Reynolds equations are then derived by substitut- 
ing these variables into the hydrodynamical equations and 
averaging over a time greater than the pulsational period. 
It is assumed that the mean flow still satisfies the hydrody- 
namical equations. Then, subtracting the mean flow equa- 
tion from the corresponding equation for mean plus fluctu- 
ating parts yi elds an equatio n for the fluctuating part. (See 
Section 1.3 of lMcComblll99Cl . which the following discussion 
parallels.) 

For example, the spherically symmetric steady state 
continuity equation is, 



d_ 

dr 



7^\ r <P U > 



0. 



(Al) 



The mean flow continuity equation is, 
l(r*<p>u)=0 



(A2) 



(see eq. (1)), and the fluctuation continuity equation is, 
dr 



■rr-ir 2 < p >U + r 2 « p> a > +r 2 < pa > I = 0. (A3) 



< p > U = 0, << p > a >=< p X a >= 0, (A4) 
so we have, 

< pa >= 0. (A5) 

The steady velocity equation is, 



du _ -1 8P_ _ d <a' r > _ GM 
dr p dr dr r 2 



(A6) 



where in contrast to equation (2) we have included 
a mean stress gradien t, or wave pressure, term (e.g., 
lLandau fc Lifshitzlll959l) . After we subtract the mean flow 
version of this equation, we get the following equation for 
the fluctuating velocity (before averaging), 



,dU 

u « — h 
dr 



or 



1 



+ 



<") r 



< a;, > 



-1 dP' 
<p> dr 



(A7) 



d<P> 



On averaging this last equation we find, 



; h hiqher order terms. (A8) 

<p> <p> or ^ *• < 

As before, all of the first order terms will vanish on aver- 
aging, leaving only thermal pressure and density terms at 
second order, which we assume are negligible. 

Equation (4) of Section 2.2.2 is the total energy equa- 
tion. Writing each variable in terms of its mean and fluc- 
tuating parts, and averaging with the assumption that odd 
order moments average to zero, yields equation (5). With 
the additional approximations described in that section, this 
equation reduces to equation (6). 

The various approximations provide a simple closure 
of the Reynolds moment equations in this application, but 
they do so via the assumption that dissipation in the wind 
formation region has taken the system to a relaxed state 
characterized by small fluctuations. Therefore, the equations 
do not provide a valid description of the unrelaxed dynamics 
at the onset of pulsations, nor of the strong shock zone shown 
by numerical models to exist at smaller radii (see Fig. 1-3). 



APPENDIX B: APPENDIX B: SPHERICAL 
SHOCKS PROPAGATING DOWN A POWER 
LAW DENSITY GRADIENT 

As described in the text, the numerical models show that the 
extended atmospheres and winds of LPVs relax to power law 
density and pressure profiles, and in fact, the main result of 
this paper is on how these profiles can be accounted for with 
the aid of shock pumping and mean entropy gradients. These 
conclusions are founded on the assumption that shocks prop- 
agate down the relaxed power law pressure gradient without 
significantly accelerating (or decelerating). In this appendix 
we briefly review some classical results that demonstrate the 
existence of profiles for constant shock propagation. These 
results do not seem to be well-known in the astrophysi cal lit- 
eratur e, though they are summarized in the text of Whitham 
Jl974h . 

The first result concerns the outward propagation of 
a spherically symmetric shock wave in a constant density 
medium. The wave is not necessarily a strong blast as as- 
sumed in the well-known Sedov- Taylor solution. An approx- 
imate si milarity solution to this more general cas e was pre- 
sented bv lGuderlevI J1942I) . also see Section 6.16 of lWhithaml 
il974h . Specifically, in the spherical adiabatic case the post- 
shock velocity is found to vary as, 
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pocr- u ' aub , uocr""'™. (Bl) 

Density gradients affect shocks in a manner similar to 
geometr i c conv ergence or divergence, and a related result o f 
ISakurail Jl960l) . described in Section 8.2 of lWhitharr] Jl974T) 
is useful. That is, the velocity of a planar adiabatic shock 
on a density gradient, p oc x a , varies as, 

u oc x~\ with A = a/3, (3 ~ 0.236, (B2) 

and thus, A ~ —0.472 for a — —2. Therefore, combining 
these two effects, we expect the variation of the spherical 
shock velocity on the density gradient p oc r~ 2 , to be 
roughly, 

0.472 -0.453 0.02 /T><^ 

u oc r r ~ r . (B3) 

That is, the shock velocity is not exactly constant, but it is 
essentially so within the accuracy of the approximations. 



